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Ground-state properties of a two-dimensional system of 
superconducting vortices in the presence of a periodic array 
of strong pinning centers are studied analytically and numer- 
ically. The ground states of the vortex system at different 
filling ratios are obtained using a simple geometric argument 
under the assumption that the penetration depth is much 
smaller than the spacing of the pin lattice. The results of 
this calculation are confirmed by numerical studies in which 
simulated annealing is used to locate the ground states of 
the vortex system. The zero-temperature equilibrium mag- 
netization as a function of the applied field is obtained by 
numerically calculating the energy of the ground state for a 
large number of closely spaced filling ratios. The results show 
interesting commensurability effects such as plateaus in the 
B — H diagram at simple fractional filling ratios. 
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I. INTRODUCTION 

In the mixed phase of type-II superconductors, mag- 
netic flux penetrates the sample in the form of quantized 
vortex lines^. The amount of flux carried by each vortex 
line is equal to the basic flux quantum $o — hc/(2e) — 
2.07 X 10~^G.cm^. These vortex lines form a special phys- 
ical system known as "vortex matter" . In the absence of 
any pinning sites in the material, the vortex lines form a 
triangular lattice known as the Abrikosov lattice^. 

Equilibrium and transport properties of the mixed 
phase of type-II superconductors are strongly affected by 
the presence of pinning centers, either intrinsic to the sys- 
tem or artiflcially generated. Understanding the effects 
of pinning in these systems is very important for practi- 
cal applications because the presence of pinning strongly 
influences the value of the critical current in the mixed 
phase. 

In recent years, a variety of nanofabrication techniques 
have been used to create periodic arrays of pinning cen- 
ters in thin-fllm superconductors 3~i3^ Such arrays may 
consist of microholes ("antidots")'^"'', defects produced 
by the bombardment of ion^ or electron^ beams, or mag- 
netic dots^^"^"^. These pinning centers are "strong" in 
the sense that each pinning site can trap one or more 
vortices at low temperatures. The effects of periodic pin- 
ning depend strongly on the relative values of and B, 
where Bff, = Pp$o {Pp is the areal density of the pinning 
centers) is the so-called "matching fleld", and B is the 



magnetic induction that determines the areal density po 
of vortices {po = B/^q). The flUing ratio, n, deflned as 
n = B / B^, measures the commensurability of the vortex 
system with the underlying pin lattice. The interplay be- 
tween the lattice constant of the pin array (determined 
by B^) and the intervortex separation (determined by B) 
can lead to a variety of interesting effects in such systems. 

Some of these effects have been observed in recent 
experiments. Imaging experiments using various tech- 
niques such as Lorentz microscopy* and scanning Hall 
probe microscopy^ have shown the formation of or- 
dered structures of the vortex system at low temperatures 
for commensurate values of n. Magnetization measure- 
ments'^'^'^ in the irreversible (vortex solid) regime have 
demonstrated the occurrence of anomalies at certain har- 
monics of B^. The effectiveness of pinning at integral 
values of n has been found^'^°~^'^ to produce regularly 
spaced sharp minima in the resistivity versus fleld curve. 
A pinning-induced reconflguration of the vortex lattice 
has been observed^'^ in a thin-fllm superconductor with 
a rectangular array of magnetic dots. Some of these ef- 
fects have been studied theoretically, using analytic^^ and 
numerical^^~^° methods. Experimental realizations of a 
system of interacting "particles" in the presence of an 
external periodic potential are also obtained in colloidal 
suspensions in interfering laser flelds^^, and in periodic 
arrays of optical traps^^ . 

In this paper, we have used analytic and numerical 
methods to analyze the zero-temperature structure of 
vortex arrays in the presence of periodic pinning. We 
have also carried out a numerical study of the zero- 
temperature equilibrium magnetization of a supercon- 
ducting fllm with a square array of pinning centers as 
a function of the applied fleld. In section II, we consider 
the ground states of a vortex system in a square array 
of pinning centers for flUings less than unity. We look 
at a class of structures that are Bravais lattices with one 
vortex per basis if the flUing n is of the form and 
with p vortices per basis \i n= p/q (p and q are integers 
greater than unity, with p < q). The structure with the 
lowest energy in this class can be obtained rather eas- 
ily. We flnd that the "ground-state" structure obtained 
this way matches those obtained from experiments* and 
simulations^^'^* for a large number of simple fractional 
values of n. The results obtained in this section can 
also be used to predict the ground-state structures for 
1 < n < 2. In section III, we consider the ground-state 
structures for flUings greater than two. In these calcu- 
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lations, we use simple geometric arguments to arrive at 
the ground states. This analysis is performed under the 
assumption that the range of the intervortex interaction, 
which is set by the penetration depth, is much smaller 
than the spacing between the pinning sites. We show 
that the ground-state structures obtained from this sim- 
ple analysis match the ones obtained from simulated an- 
nealing. This analysis is extended to rectangular and 
triangular pin lattices in section IV. In section V, the 
zero-temperature equilibrium magnetization of a vortex 
system in a square array of pinning centers is obtained 
by first calculating the ground-state energy as a function 
of the magnetic induction and then finding the applied 
field from a numerical differentiation of the data. The 
ground-state energies for different values of the magnetic 
induction are obtained using a simulated annealing pro- 
cedure. The calculated B — H curve exhibits interesting 
commensurability effects, manifested as plateaus occur- 
ring at simple rational values of the filling fraction n. The 
main results of our study are summarized in section VI. 

II. GROUND STATES FOR A SQUARE PIN 
ARRAY WITH FILLING RATIO LESS THAN 
ONE 

We consider a superconducting film that has a square 
array of pinning sites with lattice constant d. The mag- 
netic field is assumed to be perpendicular to the sur- 
face of the film. The "matching field" i3$ is then given 
by = <^q/cP', and the filling fraction n is given by 
n = B/B^ = B(P/^o where B is the magnetic induction. 
We assume that the pinning potential is much stronger 
than the intervortex interaction, but is of extremely short 
range. The large strength of the pinning potential implies 
that the vortices must occupy pinning sites as long as the 
number of vortices does not exceed the number of pinning 
sites. We also assume that a pinning site can not accom- 
modate more than one vortex. If the pinning centers in 
the film are microholes, then this assumption amounts 
to the requirement that the radius of each hole is close 
to two times the coherence length ^ of the superconduc- 
tor. These assumptions ensure that interstitial vortices 
appear only when the filling fraction n is greater than 
unity. The assumed short range of the pinning potential 
can be justified if the defect diameter is small compared 
to the defect spacing d. Another assumption that we will 
make in most of our calculations is that the intervortex 
interaction falls off rapidly with distance. This is ensured 
if the penetration depth A is much smaller that the pin- 
lattice spacing d. In our calculations, we take the ratio 
A/d to be 10. This value is appropriate for the pin lattice 
of Ref. 8. We consider temperatures that are low enough 
to neglect effects of depinning and vortex-lattice melting. 
The problem of finding the structure of the vortex system 
then reduces to locating the ground state in the presence 
of the pinning potential. 



Consider now fillings of the form n = 1/q, q being 
an integer greater than one. Let us look at Bravais lat- 
tices that can be formed for a specific n by distributing 
the vortices on the square pin lattice with one vortex 
per basis. The motivation for considering such lattices is 
that this will automatically ensure that there is no shear 
of the vortices with respect to the pin lattice, since the 
forces on a vortex due to other vortices will add up ex- 
actly to zero. The unit cell area of these structures has 
to be q.d^. So the possible unit cells can be obtained by 
factorizing q into products of the form r.s (r and s are 
integers) , arranging the vortices at the corners of rectan- 
gles of dimension rd x sd, and then sliding the parallel 
sides relative to each other. This procedure produces a 
large number of structures depending on the value of n 
and we have to pick the one that minimizes the energy. 
For small values of q, this can be done by hand, but as 
q becomes large and highly factorizable, the number of 
possible structures increases rapidly. In such cases, we 
have resorted to computers to generate these structures 
and compare their energies. 

The structures so obtained for fillings 1/2 and 1/4 
match those found in the imaging experiment^. Also 
for fillings 1/2, 1/3, 1/4, 1/5, 1/8, 1/9, 1/10 and 1/15, we 
find the same structures as those obtained by solving 
the "greedy lattice gas model" exactly. This is un- 
derstandable because when the intervortex interaction 
falls off rapidly as the distance is increased above the 
defect spacing d, the ground state can be attained by 
finding the lattice that maximizes the shortest distance 
between vortex pairs. If two structures have same value 
and number of shortest distances, then the next shortest 
distance should be maximized, and so on. For fillings 
1/2, 1/3, 1/4, and 1/5, our analysis also yields the same 
structures as those found in the large Ue [U,, is the en- 
ergy of on-site repulsion between two electrons) limit of 
the neutral Falicov-Kimball modeP**. In Fig. l(a,b,c), 
we have shown the structures so obtained for a few fill- 
ings of the form n = l/g. The ground-state structure 
shown in Fig. 1(b) for n = 1/5 is different from that 
found in Ref. 16 from a simulated annealing calculation. 
This difference is probably due to the use of a different 
(logarithmic) intervortex potential in Ref. 16. 

Let us now compare the energies of the nearest and 
next-nearest neighbors in one of these lattices. The in- 
teraction energy between two vortices separated by a dis- 
tance r is given by the expression, 

where Kq is the zeroth order Hankel function of imagi- 
nary argument and t is the film thickness. For n = 1/2, 
the nearest neighbor distance is \/2d and the next-nearest 
neighbor distance is 2d. So the interaction energies are, 
to within a constant prefactor, given by 

UnOcKo =0.2x 10-^ 
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U„n cx Ko l^ — j = 0.6 X 10-1°. 

One can see here that there is orders of magnitude dif- 
ference in these energies which cannot be compensated 
by differences arising from interactions with more distant 
neighbors. This difference is going to be more prominent 
at lower densities. This tells us that the maximization of 
the shortest intervortex distance in a lattice for a given 
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FIG. 1. The ground state structures for a few filling frac- 
tions n < 1. The different filling fractions are (a) 1/2, (b) 
1/5, (c) 1/9, (d) 2/3, (e) 3/5, and (f) 2/7. The dots in the 
figures represent the pinning sites and the circles represent 
the vortices. 

filling would lead to the ground states, provided the lat- 
tice spacing is large compared to the penetration depth 
of the film. This, in fact, is exactly the definition of the 
greedy lattice gas. However, one has to be cautious about 
this method because, as noted in Ref. 23, the structures 
can be strongly dependent on the form of the potential 
in certain ranges of n and we can even have aperiodic 
structures as ground states. The ground-state structures 
shown here have been cross checked with simulations to 
ensure that they are indeed the lowest energy configura- 
tions. To give an example of a case where this treatment 
does not lead to the true ground state, we found that for 
filling 1/16, the energy per vortex for the structure with 
minimum energy obtained this way was greater than that 
for filling 1/15, implying that the structure obtained for 
n= 1/16 was not the ground state. 

When the filling fraction is of the form p/q with p not 
equal to one, one can look for ground states in a subset of 
structures where the unit cell has size qcP with p vortices 



in a basis. We have shown in Fig. l(d,c,f) some of the 
ground-state structures obtained this way. These struc- 
tures match those obtained from our simulated anneal- 
ing calculation. These ground states show the "stripe" 
structure predicted by Watson^^ and Kennedy^^ in ap- 
propriate density ranges. 



III. GROUND STATES FOR A SQUARE PIN 
ARRAY WITH FILLING RATIO GREATER 
THAN TWO 

If n is greater that 1 but smaller than 2, then the 
ground-state structures are similar to the ones for the 
case of n less than one. The only difference is that 
the pinning sites are all occupied and the centers of the 
square unit cells of the pin lattice now act as new pinning 
centers for interstitial vortices. But things look different 
when the filling goes above the value n = 2. For such val- 
ues of n, wc can no longer place the interstitial vortices at 
the centers of the squares and look for simple structures 
obtained this way. Also, we now have to start looking 
into the stability of the structures since the square sym- 
metry would not be present. 

A. The ground state for n = 5/2 

Here we arc faced with the task of placing more than 
one vortex in a square. Before going to the problem of 
finding the ground state for n = 5/2, let us ask a more 
basic question: given a single unit cell of the square pin 
lattice with each corner occupied by a vortex, how can 
we arrange two more vortices inside this square so as 
to minimize the energy. Since "greedy lattice gas" has 
been a good approximation for the previous cases, we try 
to tackle this problem by "maximizing the shortest dis- 
tance" method. In order to stabilize an interstitial vortex 
by maximizing the shortest distance, its distance from at 
least three nearest vortices must be the shortest distance. 
It is also required that these nearest set of vortices must 
be spread in such a way that if wc draw straight lines 
from the vortex in question to these neighbors, the an- 
gles formed by adjacent lines must be less than 180°. The 
proof of this statement is given in Appendix A. 

It can be seen from the symmetry of the problem that 
we have to place the two vortices on the lines joining 
the centers of the sides to meet the condition mentioned 
above. This leaves us with only two possible ways of 
doing it, which are shown in Fig. 2. In (a) the shortest 
distance D^a can be obtained by solving the equation 

i(d-Z)^,f + ^ = z?L- (2) 

On solving this equation, we get AB = BC = BE = 
Dma = (\/7- l)d/3. In (b) the vortices A',B' and C" 
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FIG. 2. Putting two vortices in a square. The configuration 
in (a) is the global energy minimum, and the configuration 
in (b) can at best be a local minimum of the energy. The 
angles 6 — 15° and 4> — 24.49°. The square is drawn for easy 
visualization and the dotted lines are the bisectors of the sides. 
Note that there is already one pinned vortex at each of the 
corners of the square. The distances AB = BC = BE = Dma 
and A'B' = B'C = B'E' = Dmb- 



form an equilateral triangle. Thus the nearest neighbor 
distance A'B' = B'C = B'E' is 



£>„6 = sec(15°)| 



(3) 



The angle (f> in (a) is 24.49°, and the angle 9 in (b) is 
15°. The interaction energies corresponding to these two 
distances for d/A = 10 are 

Uab oc Ko t^^^^d] ~ 2.2 X 10"^ 



Ua'b' oc Ko 



sec(15°) 
2A 



3.0 X 10" 



From comparing these two energies it is clear that (a) 
is the global minimum, whereas configuration (b) can at 
best be a local minimum. 

Coming back to the n = 5/2 case, we now have to build 
up the lattice with equal number of two types of squares 
- one with two interstitial vortices and the other with 
one interstitial vortex. Note that here we have neglected 
structures that have three or more interstitial vortices 
inside a square unit cell because such structures would 
drastically bring down the nearest neighbor distance. Let 
us now look at the possible units cells of size 2d x 2d that 
can be made out of these two types of squares. These are 
shown in Fig. 3. If one constructs the lattice with these 
unit cells, configuration (b) offers the least number of 
next-nearest neighbors, the number of nearest neighbors 
being the same in all the cases. So one can expect (b) to 
be the ground state unit cell, at least for large values of 
d/X. However, unit cell (a) is preferred if d/X is not very 
large. This can be understood in the following way: the 
advantage that (b) has over (a) is that it has only half the 
number of next-nearest neighbors (interactions like that 




FIG. 3. Possible 2d x 2d unit cells for n = 5/2. The unit 
cell (a) is the lowest energy configuration for d/X = 10. For 
much larger values of d/A, the unit cell with the lowest energy 
will be the one in (b). 



between vortices m and n) compared to (a). But this 
is done at the cost of bringing in interactions like those 
between vortex pairs (p, q) and (p, t) for every "gain" of a 
next-nearest neighbor interaction. The energies of these 
two interactions for d/X = 10 are found to be quite close. 
These energies are 



Umn OC Ko{rmn/X) ~ 3.2 X 10" 



Upq oc Ko{rpq/X) ~ 1.9 X 10" 



where 



and 



are the distances between vortices m 



and n, and p and q in Fig. 3, respectively. It is clear from 
this comparison that the unit cell (a) would be preferred 
for d/X = 10. 

The ratio of interaction energies of the next-nearest 
and the nearest neighbors is 0.07 for this lattice when 
d/X = 10. This energy difference is appreciable here, so 
that we can expect that the unit cell we arrived at is 
the correct one. Note that any net force that might be 
present on one of the interstitial vortices due to the asym- 
metry in the structure can be compensated by extremely 
small displacements from the positions obtained from the 
"maximization of the shortest distance" method. 



B. The ground state for n = 3 

When the filling fraction equals 3, we have to build 
up the lattice using blocks of type (a) in Fig 2. Again, 
looking at unit cells of size 2d x 2d or smaller, we have 
the configurations shown in Fig. 4 to consider. Here it is 
easy to see that unit cell (b) is preferred over the others. 
This is because it is the configuration that maximizes the 
minimum distance between any two vortices in different 
squares, the distances between vortices within one square 
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FIG. 4. Possible 2d x 2d unit cells for n — 3. Configura- 
tion (b) can be easily seen to offer the maximum next-nearest 
neighbor distance, the nearest neighbor distances being the 
same and equal in number in each case. 



being the same in all the configurations. Again compar- 
ing the nearest interaction and the next-nearest one, we 
have 



Un oc Ko {vsp/X) 2.2 X 10" 



Unn OC Ko {rpJX) ~ 1.9 X 10-^ 

There is appreciable difference between these two val- 
ues, and hence, the ground state we have arrived at is 
reasonable. When the filling lies between 2 and 3, one 
can safely assume that the ground state structure can be 
built up using squares of type (a) in Fig. 2, and squares 
that have one vortex at the center. In fact we make use 
of this in our simulations to arrive at the ground states, 
as described in section V. 



C. The ground state for n = 4 

Here we have to place three interstitial vortices in one 
square. This is nontrivial since even if we ensure that the 
shortest distance is maximized in one sqiiarc, two vortices 
in nearby squares may be closer to each other than the 
shortest distance within a square when we create the lat- 
tice. Note that we did not come across this problem in 
the n = 5/2orn = 3 fillings. To illustrate this problem, 
we show in Fig. 5(a) a single cell of a square lattice with 
three interstitial vortices, which may very well be a lo- 
cal minimum configuration. But if we try to build the 
lattice using this cell, we can not do it without bringing 
the vortices in nearby squares closer than the minimum 
distance in an individual cell. So what we need to do is 
to look for a pattern that will include vortices in differ- 
ent squares while doing the minimization of the shortest 
distance. We can solve this problem trigonometrically. 
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FIG. 5. Two possible ways of arranging three vortices in a 
square so as to maximize the shortest distance. The configu- 
ration in (a) offers the best arrangement if looked in isolation, 

but the configuration in (b) wins out when one has to con- 
struct a lattice of the unit cells. The angles = 15° and 
ip = 60". 



Consider the figure shown in Fig. 6. The solution we are 
looking for can be obtained by solving the equations 



AB = BC = CA = D, 



CB' = AC" = AP = D,. 



(4) 



(5) 



Note that here we have assumed a unit cell size d x d. 
On solving these equations, we obtain the unit cell shown 
in Fig. 5(b). In the figure, the angle 6 equals 15° and 
the angle (p is 60°. In this lattice, the shortest distance 
is Ds = sec(15°)(i/2 and the next shortest distance is 
Dns = 3dsec(15°)/(4\/2). 

This simple solution may not be the correct one if the 
lattice spacing is not large enough. For example, if d/A = 
10, as we have been assuming when comparing energies, 
then the nearest and next-nearest neighbor interaction 
energies turn out to be really close. Hence we can not 
rule out the possibility of the lattice arranging in such a 
way that the shortest distance is reduced so as to decrease 
the number of nearest or next-nearest neighbors. The 
relevant energies for d/X = 10 are 

Ko (^^^ ~ 3.1 X 10-^ 

Ko (^^^ ~ 2.1 X 10-3. 

So the ground state obtained above is guaranteed to be 
the correct one only for much larger values of d/X. The 
ground states we have obtained for n = 5/2 and n = 3 
match with the images from experiments* and also the 
results of simulations-'^^. But for n = 4, the structures 
found in experiments and simulations are different from 
the one shown in Figs 5(b) and 6. This is expected, since 
in the experiment the value of the ratio d/X was close 
to ten. Our simulation with d/X= 10 gives the ground- 
state structure show in Fig. 7 (a), which is similar to the 
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the longer side of the rectangle. We shall consider here 
only the cases where the filling is greater than one. In 
the absence of the square symmetry, it is obvious that 
the ground state structure will depend not only on the 
penetration depth A, but also on the ratio l/h. In the fol- 
lowing analysis, we shall always assume that A is much 
smaller than b, the shorter side of the basic rectangular 
pinning cell. When the filling is 2, for values of l/h less 
than v^. the ground state is one where each interstitial 
vortex is at the center of the rectangle, since this ensures 
that the shortest distance is maximized (see Fig. 8(b)). 
But when the aspect ratio exceeds \/3 and the interstitial 



FIG. 6. The distances one has to equate when maximizing 
the shortest distance for n = 4, taking into consideration 
distances between vortices in neighboring squares. Note that 
the unit cell size here is one square unit. 



one obtained in the experiment^. In our simulations with 
very large values of d/\, we get structures similar to that 
in Fig. 6. The simulation result for rf/A = 50 is shown 
in Fig. 7 (b), which matches well with the predicted 
structure. It is to be noted that the simulation result 
was obtained by starting the system near the expected 
ground state. So the claim is that it offers at least a 
local minimum of the interaction energy. The simulations 
were carried out for different system sizes form 2d x 2d 
to lOd X lOd to rule out any dependence of the results on 
the boundary condition. 



30<> 
20" 
10<> 



10 20 30 

(a) 



40 



200' 
150' 
100<t' 
50> 
04 



50 100 150 200 

(b) 



FIG. 7. The ground state for n = 4, obtained by simu- 
lated annealing when the ratio of the penetration depth to 
the pin-lattice spacing is (a) 10 and (b) 50. The dark dots 
denote the pinning sites and the axis labels are in units of the 
penetration depth. 



IV. GROUND STATES FOR RECTANGULAR 
AND TRIANGULAR PIN ARRAYS 

One can extend this type of analysis to pinning arrays 
with other symmetries for finding the least energy struc- 
tures for simple filling fractions. Let us first consider the 
case of a rectangular array of pinning sites with a pin- 
ning unit cell of dimensions I x 6, where we take I to be 
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(b) 



FIG. 8. The unit cells of the ground state, obtained 
by maximizing the shortest distance, when the filling is 
2 for a rectangular array of pinning sites, (a) The unit 
cell when the aspect ratio is greater that Vs. The inter- 
stitial vortex is displaced horizontally from the center of 
the rectangle by distance Ds- In the figure, the distances 
AB = AC = AD = BF = BE. (b) The ground state unit 
cell when the aspect ratio is less that \/3, where the intersti- 
tial vortex is at the center of the rectangle. 

vortices are placed at the centers of the rectangles, the 
distance between two interstitial vortices in neighboring 
cells would be shorter than that between an interstitial 
vortex and the closest pinned vortex. This would lead to 
a displacement of the interstitial vortices sideways from 
the center, along the bisector of the shorter sides of the 
rectangle, to maximize the shortest distance. The result- 
ing structure is shown in Fig. 8 (a). The displacement 
of the vortex from the center is given by 



D. = 



-I + - 962 
6 ■ 



(6) 



It is worth noting that since the vortex in the center 
would be moving towards two of the pinned vortices, and 
away from only one interstitial vortex per unit cell, the 
displacement will approach the value given above only 
when the ratio l/b is appreciably large and in the limit 
of small penetration depth compared to the sides of the 
rectangle. For example we have found in our simulations 
that even when the ratio l/b is 2, the ground state for 
6/ A = 15 is one where the interstitial vortex is very close 
to the center, whereas for the l/b = 3 and b/X = 15, the 
ground-state structure is quite close to the one obtained 
from maximizing the shortest distance. Some of our sim- 
ulation results for filling equal to two are shown in Fig. 
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12 (a,b). Also, if the ratio l/h becomes too large, the 
analysis will have to include more than two of the inter- 
stitial vortices, since now the solution like that shown in 
Fig. 8 (a) can lead to two vortices being closer in the 
next nearest cells or ones even further apart. 

In trying to arrive at the lowest energy structures for 
fillings 5/2 and 3, it is important to determine how one 
can accommodate two vortices in a rectangular cell with 
the shortest distance being maximized. There are two 
possible minima that one has to consider: one where the 
vortices are arranged along the line dividing the shorter 
sides, and one where they are arranged along the line 
dividing the longer sides, as shown in Figs 9 (a) and 9 
(b), respectively. The shortest distance in each case is 
given by, 



D 



V4P + 2,lP-l 

3 ' 
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If one considers the distances within the cell, configura- 




(b) 



FIG. 9. The two possible ways by which the shortest dis- 
tance can be maximized when there are two vortices in a 
rectangular cell, (a) In this case the vortices arc placed on 
the line that divides the shorter sides of the rectangle. The 
distances AB = AC = AD = DE = Dsi. (b) Here the vor- 
tices are placed parallel to the shorter sides and the distances 
A'B = A'E = A'D' = CD' = Ds2. For l/b less than 2, the 
configuration in (b) leads to a larger shortest distance (within 
the cell) than the one in (a). If Z/& is greater than 2, then the 
distance D32 becomes greater than b and the vortices spill 
over into the next cell. 

tion (b) in Fig. 9 gives the lowest energy. But for large 
values of the ratio l/b, this configuration is disfavored 
since it allows the vortices in one rectangle to get close 
to those in a neighboring one. Also, for l/b > 2, the 
interstitial vortices "spill over" into the next cell, since 
the distance Ds2 becomes greater than b. So, one has 
to work out the structures for hllings 5/2 and 3 case by 
case. For fillings between 2 and 3, one has to choose 
appropriate number of two- vortex rectangles of the right 
kind and single-vortex ones and arrange them so as to 
maximize the shortest distance appearing in the struc- 
ture. We have looked at the 2x2 unit cells possible for 
filling 5/2 for two values of the aspect ratio, l/b = A and 



l/b = 5/4. The unit cells that provide the largest min- 
imum distance are shown in Fig. 10. Note that when 
l/b = 4, the vortices are arranged parallel to the longer 
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FIG. 10. The unit cells for n — 5/2 for a rectangu- 
lar array of pinning sites, (a) The unit cell when the as- 
pect ratio is 4. The distances AB = AC = AE = BD. 
(b) The unit cell when the aspect ratio is 5/4. Here 
PV = PU = PQ = QW = QS = D,2 and QR= RS = RT. 
When this structure is repeated periodically, the "image" of 
the vortex at P would be at the same distance QR from R. 
This would ensure the stability of the vortex at -R. 

side and in the other case, parallel to the shorter side. 
Also when l/b = 5/4, the vortex in the single- vortex 
rectangular cell is not located at the center, but slightly 
displaced sideways along the bisector of the shorter sides 
to facilitate the maximization of the second shortest dis- 
tance involved. One should again keep in mind that this 
sort of analysis will go wrong if the aspect ratio is too 
large, since then the distances between vortices in next- 
nearest or further neighbor cells will become important. 
In Fig. 12 (c) we have shown the ground state struc- 
ture obtained from simulations for aspect ratio 4, with 
b/X = 5 and in Fig. 12 (d), the ground state for aspect 
ratio 5/4, with b/X = 24. In Fig. 12 (d), one can notice 
the slight shift of the vortex from the center in rectangu- 
lar units having a single interstitial vortex. The shift is 
not as large as expected from the above analysis. This is 
due to the large but finite value oi b/X. 

For filling equal to 3, the lowest energy structures ob- 
tained by considering 2x2 unit cells for two values of 
the aspect ratio, 2 and 5/4, are shown in Fig. 11. Here 
too, for large values of the aspect ratio, the structure is 
composed of rectangular cells where the interstitial vor- 
tices are aligned parallel to the longer side (Fig. 11(a)), 
whereas when the aspect ratio is smaller, the structure 
is made up of an alternating arrangement of rectangular 
cells of both types (Fig. 11(b)). The simulated anneahng 
results for similar values of the aspect ratio (Fig. 12 (e,f)) 
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FIG. 11. Unit cells for n = 3 for a rectangular array of 
pins, (a) The unit cell when the aspect ratio is 2. The 
distances AB = AC = AE = BD = BF. (b) The unit 
cell for the same filling but for aspect ratio 5/4. Here the 
distances PQ = PS = PT = QR = QU = Ds2 and 
TV = UV = VW = WY = WX = Dsi. 

yield the structures obtained from the above analysis. 

Ground state structures obtained by simulated anneal- 
ing for a rectangular pin array with l/b = 2 and integral 
values of n are reported in Ref. 17. In that study, the 
intervortex interaction was assumed to depend logarith- 
mically on the intervortex distance. The ground state 
structure found in Ref. 17 for n = 2 is similar to that 
shown in Fig. 12(b), but the structure found there for 
n = 3 is quite different from those of Fig. 12(e,f). This is 
another example of the importance of the detailed nature 
of the intervortex interaction in determining the struc- 
ture of the ground state. 

For a triangular array of pinning sites, it is easy to see 
that when the filling is greater than 1 and less than 3, 
the interstitial vortices will be placed on the centroids of 
the triangles in the limit where one can safely apply the 
method of maximization of the shortest distance. So the 
ground states when n is between 1 and 3 will be made 
up of parallelogram cells of the form shown in Fig. 13. 
These unit cells match well with the results of molecular 
dynamics simulations^^'^^. 



V. EQUILIBRIUM MAGNETIZATION OF THE 
GROUND STATES 

In this section, we describe a calculation of the zero- 
temperature, equilibrium magnetization of a thin- film su- 
perconductor in the presence of a square array of pinning 
sites. The region in the B — H plain we are interested in 
is that just above Hci, when the flux tubes start enter- 
ing the sample. The idea is to find the free energy F of 
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FIG. 12. The ground state structures obtained by simu- 
lated annealing when the pinning array has rectangular sym- 
metry. The parameters for different plots are, (a) n = 2, 
l/b = 3 and b/\ = 15, (b) n = 2, l/b = 4/3 and b/\ = 15, 
(c) n = 2.5, l/b = 4 and b/X = 5, (d) n = 2.5, l/b = 5/4 and 
b/X = 24, (e) n = 3, l/b = 2 and b/X = 10 and (f) n = 3, 
l/b — 5/4 and b/X = 16. The dots denote the pinning centers 
and the circles represent the vortices. Note that the box sizes 
are not to scale. 




FIG. 13. Basic building blocks for generating ground states 
for triangular pinning arrays when the filling is between 1 and 
3. (a) When there is a single interstitial vortex in a parallel- 
ogram, the shortest distance can be maximized by placing it 
at one of the two centroids of the triangles involved. The dis- 
tances AB — AD = AC = a/^/3, where a is the length of the 
side of a pinning cell, (b) When two vortices are to be placed 
in a single pin cell, they have to be at the two centroids. 
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the ground state as a function of the magnetic induction 
-B, and then to obtain the apphed magnetic field H by 
taking a derivative of the free energy with respect to B. 
Since we are considering the zero-temperature case, the 
free energy is just the internal energy of the flux lattice. 
Since we are looking for a nearly continuous variation of 
the internal energy for taking the derivative, we need to 
locate the ground states for filling fractions separated by 
small intervals. This would be difficult to do analytically, 
since the unit cells for some filling fractions can be arbi- 
trarily large. Also, as n becomes large, the simple proce- 
dure of maximization of the shortest distance is not going 
to yield the correct ground-state structures. So we have 
resorted to simulations to determine the ground states. 
In particular, we have used the simulated annealing tech- 
nique to locate the global minima (or at least low-lying 
local minima close in energy to the global ones) of the 
part of the internal energy associated with intervortex 
interactions. 

The Hclmholtz free energy per unit volume of the su- 
perconductor at zero temperature in the presence of the 
pinning sites is 



Fs{n) 



En 



<P 



(9) 



where the first term is the line energy, the second term is 
the interaction energy and the third term is the pinning 
energy. Here, ei is the line energy per unit length, Cp is the 
pinning energy per unit length and En is the interaction 
energy per unit volume for filling fraction n. We note 
here that the pinning energy increases linearly with n till 
n becomes one and then remains constant, since multiple 
occupation of a pinning center is not allowed. Further, 
for simplicity, we express the pinning energy as 



gp = me; 



(10) 



where m is a positive number whose magnitude depends 
on the nature of pinning. The interaction part of the free 
energy. En, is the computational input. Once we know 
the free energy, we can compute the applied magnetic 
field H using the relation 
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Using the standard expression-'^ for e; and taking the log- 
arithm of the Ginzburg-Landau number to be equal to 2, 
we get the following expression for the applied field as a 
function of the filling fraction: 



H = ,^[l-m{l- Q{n -!)}]+ " 
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Here E'n is given by the expression 
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FIG. 14. The ground states for (a) n = 9/4, (b) n = 5/2, 
(c) n = 11/4, (d) n = 3 obtained by simulated annealing, 
as discussed in the text. The unit-cell size is 4d x Ad. The 
dark dots denote the pinning sites and the circles denote the 
vortices. The axis labels are in units of the penetration depth. 



where N is the number of basic pinning squares in the 

system and is the separation between vortices i and j 
in the ground state for the filling fraction n. 

The size of the systems we simulated varies from 2dx2d 
to 8d X 8d. In all cases, we used periodic boundary con- 
ditions to minimize surface effects. So the minimum 
difference between two consecutive filling fractions was 
An = 1/64. The ratio d/X was taken to be 10, as 
in our previous analysis. In order to save computation 
time, the vortices were allowed to stay only at the pin- 
ning sites when the filling was less than one. For fill- 
ings between one and two, every pinning site was occu- 
pied by a vortex which was never moved and the extra 
ones were allowed to move near the centers of the ba- 
sic pinning squares. When the filling was greater that 
two and less than three, the vortex configurations were 
constructed using basic units of squares containing one 
vortex at its center and squares containing two vortices 
placed such that the shortest distance within a square 
is maximized (as in Fig. 2(a)). These units were then 
moved around and twisted while cooling to arrive at the 
minimum energy states. This procedure helped us to 
track low-lying minima faster than if we allowed vortices 
to move freely. Once the basic striicture was thus ob- 
tained, the vortices were allowed to move freely during a 
second cooling schedule starting from a lower tempera- 
ture to obtain the lowest-energy structure. In Fig. 14 we 
have shown some of the ground state structures we have 
obtained this way for fillings between 2 and 3. For fill- 
ings 5/2 and 3, we find that the structures match those 
obtained in experiments^, as well as in our analysis using 
maximization of the shortest distance. The structures 
for n = 9/4 and n = 11/4 may not be the actual ground 
states, either due to the smallness of the unit cell of our 
simulation or due to the presence of many nearly degen- 
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FIG. 15. The energy of the vortex lattice in the presence 

of a square array of pinning sites (upper dashed curve) . The 
energy of the triangular vortex lattice in the absence of pin- 
ning is also shown for comparison (lower solid curve). The 
system size is 8c? x 8d. The abscissa is the filling fraction n 
and the ordinate is the total energy per unit thickness in units 
of ^^/(STr^A^). 

In Fig. 15 we have plotted the ground state energies 
obtained from the simulation for different fillings. The 
simulation unit cell was 8d x 8d and the energies were 
computed for fillings 1/64 to 3. The upper curve shows 
the results ol)tained in the presence of the pinning sites 
and the lower curve is the energy of the triangular lattice 
for the same density of vortices. Note that we have not 
included the pinning energy in the plot. This would bring 
down the upper curve below the curve for the pin-free 
case. 



3 
2.5 

2 



§ 1.5 



1 

0.5 



15.75 15.85 15.95 

H/B, 



16.05 



FIG. 16. The dependence of the magnetic induction B on 
the applied field H in the entire region of our simulation. A 
number of plateaus can be seen at points corresponding to 
simple rational filling fractions mentioned in the text. Only 
the vertical lines correspond to the data points obtained from 
our calculation; the dotted lines are guides to the eye. Both 
B and H have been scaled by the matching field B^. 



From the energy versus filling fraction data, one can 
find the applied field using Eq.(ll) and then compute the 
magnetization M using the relation 



B = H + A'kM 



(14) 



In Fig. 16 we have plotted B versus H in the entire 
range of filling for which simulations were carried out, 
from n = to n = 3. Figs 17 and 18 show magnified 
versions of this plot in the regions between filling frac- 
tions and 1, and between 2 and 3, respectively. Note 
that we have not explicitly included the pinning energy 
term in our analysis. This term would just add a con- 
stant contribution to H for fillings up to 1 . The features 
of the curve from n = 1 to n = 2 are the same as those 
in the interval between n = to n = 1. This is due to 
the fact that the ground state structures are similar in 
the two regions (see section III). The B — H plot shows 
flat regions at values of B corresponding to fillings 1/8, 
1/5, 1/4, 1/2, 3/4, 4/5, 7/8, 1, 9/8, 6/5, 5/4, 3/2, 7/4, 
9/5, 2 in the filling fraction range between and 2. Also, 
in the range of n between 2 and 3, there are roughly two 
plateaus, appearing near n = 2.3 and n = 2.6. 
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FIG. 17. Expanded view of the plot of Fig. 16 in the region 
between n = B/B^ = and n = 1. One can see plateaus 
appearing at n = 1/4, n = 1/2 and n = 3/4. The dark dots 
are the data points. The light dotted line is drawn as a guide 
to the eye. 

The observed values of the filling fractions between 
n = and n = 1 at which the plateaus occur indicate 
that these values of n correspond to fillings where the 
introduction of a new vortex into the system leads to the 
appearance of a shorter distance than those existing in 
the lattice or in its dual (i.e. the lattice obtained by re- 
placing particles by holes, and vice versa) This makes 
sense because the introduction of this shorter distance 
brings in a larger energy scale, leading to a discontinu- 
ous change in the derivative of the energy with respect 
to the filling n. In the simulations, we have not scanned 
very small intervals of n. Also, for some of the fillings, 
the ground state may not have been obtained in our sim- 
ulated annealing calculation. For these reasons, we can 
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not say anything definite about the true nature of the 
B — H curve. It is possible that this curve has plateaus 
and discontinuities occurring at all scales (e.g. at all ra- 
tional values of n). The noisy nature of the B — H plot 
in the range where n lies between 2 and 3 (see Fig. 18) 
is also due to these difficulties. However, this plot shows 
clear signatures of two plateaus that appear near n = 2.3 
and n = 2.6. These can be understood as happening 
when first there is an occurrence of squares containing 
two vortices coining next to each other diagonally (as in 
Fig 3(a)), and again when they have to be next to each 
other with a common side (as in Figs 3(b) and 3(c)), as 
the value of n increases from 2 to 3. 

3| ■ 
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FIG. 18. Expanded view of the plot of Fig. 16 in the region 
between n = 2 and n = 3. One can see plateaus at fillings 
n ~ 2.3 and n ~ 2.6, as indicated by the arrows. The dark 
dots arc the actual data points and the light dotted line is 
shown as a guide to the eye. 



VI. SUMMARY AND DISCUSSIONS 

In this paper, we have reported the results of an an- 
alytic calculation of the lowest-energy states of a vortex 
system in the presence of a periodic array of strong pin- 
ning centers, each of which can trap only one vortex. 
We have considered several different lattice structures of 
the pin array and a large number of filling fractions in 
the range between zero and four. The analytic calcula- 
tions are based on the principle of maximization of the 
shortest intervortex distance. We have argued that this 
principle leads to the exact ground states when the spac- 
ing of the defect lattice is large compared to the range of 
the intervortex interaction set by the value of the pene- 
tration depth. This principle has been used, in combina- 
tion with simple geometric considerations, to obtain the 
ground states for several values of the filling fraction n. 
The ground-state structures so obtained are found to be 
identical to those found in imaging experiments® and in 
earlier simulations^^'^®. We have also carried out simu- 
lated annealing calculations of the ground states in order 



to test some of the predictions of the analytic approach. 
In all cases, we foimcl that the analytic results agree with 
those of our numerical calculation. 

We have also described the results of a numerical calcu- 
lation of the equilibrium magnetic induction B and mag- 
netization M of a planar superconductor with a square 
array of pinning centers as functions of the externally 
applied field H. We show that the interplay between 
the lattice spacing of the pin array and the intervor- 
tex separation set by the value of B leads to interest- 
ing commensurability effects, showing up as plateaus and 
discontinuities in the B vs H plot at simple rational 
values of the filling fraction n. Anomalies in the irre- 
versible magnetization of thin-film superconductors with 
periodic arrays of pinning centers have been observed at 
certain integral values of n in experiments^'^ and simu- 
lations^^. The presence of a periodic array of pins is also 
expected^^ to produce anomalies in the equilibrium mag- 
netization of the high-temperature vortex liquid at small 
integral values of n. Our results show that these com- 
mensurability effects are more pronounced in the field- 
dependence of the equilibrium magnetization and mag- 
netic induction of such systems in the low-temperature 
vortex-solid regime. Experimental investigations of these 
effects would be most welcome. 



APPENDIX A: CONDITION FOR MAXIMIZING 
THE SHORTEST DISTANCE 

Given a particle fixed at some point in a plane, the 
problem is how to place three other particles around the 
first one in such a way that if we try to move the first 
particle from its position, it will get closer to at least 
one of the three particles. The solution is as follows. If 
we draw straight lines from the particle that we want 
to move to the other particles, then each angle between 
adjacent lines must be less than 180°. In other words, 
it should not be possible to draw a straight line through 
the particle in question in such a way that all the other 
three particles lie on one side of the line. 

Let the particle that we want to move be at Pq (see Fig. 
19). Let us place two particles at Pi and P2, anywhere 
on the plane. This can be done since one of the angles 
between any two lines will always be less than or equal 
to 180°. In Fig. 19, AE and CF are the tangents to 
the circles centered at Pi and P2 and passing through 
the point Pq. The presence of these particles restricts to 
within the angle ^ the direction in which the particle at 
Pq can move without decreasing its distance to Pi and 
P2. Now let us check where we can fix the third particle 
such that the aforementioned condition is met, that is the 
particle at Pq cannot be moved without bringing it closer 
to one of the particles at Pi , P2 and P3 . We will consider 
two possible cases - the third particle in the region BPqD 
of the plane and outside this region. 

When the third particle is in the region BPqD (point 
P3 in the figure), it is not possible for the particle Pq to 
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FIG. 19. Geometry of the problem of maximizing the short- 
est distance (see text). 



move from its position without decreasing the distance to 
any one of the particles. This can be seen if one drops the 
perpendiculars to the lines AE and CF from the point 
P3. If the point Pq is moved along either of these lines, 
it will be getting closer to the point P3 . But that in turn 
implies that it cannot be moved into the region FPqE at 
all without decreasing the distance from the point P3. 

If the particle is outside the region BPqD, for exam- 
ple a point like P3, then it is easy to see that by moving 
along one of the tangents to the circles at the point Pq, 
the particle at Pq can move away from all the three other 
points. Thus we find that the earlier statement we made 
is proved. This gives us a nice way to maximize the short- 
est distance since all one has to do is to make the three 
distances involved equal, so that if the central particles 
tries to move, then one of the distances has to decrease. 
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